
% this function works for tracking vasodilation with a intravascular tracer
% like dextran-texas red
%
% this version uses kymoTrackUsingGradient to get initial edge position
% estimated, then does a refining step fitting the actual edge using
% fitedge (from the guess generated by kymoTrackUsingGradient)
%
% update from v3_0 - smoothing done with LoG instead of median to help
% emphasize edges
%
% update from v3_1 - smoothing kernel is adjusted to the current zoom
% instead of arbitrary pixels (v3_2 is for tracking just the vessel
% edges...)

function results = tracerTrackRecipe_v3_3(mov,vesselCenterIm,autoSelectCenters,pix2mic)

%% set parameter values
medFiltRadius = 3;

%% parameters for LoG filtering
kernsz = 10;     %in microns
logsig = 1.5;     %in microns
kernsz = ceil(kernsz/pix2mic);
logsig = logsig/pix2mic;
[X,Y] = meshgrid(-kernsz:kernsz, -kernsz:kernsz);
kern = (1/(pi*logsig^4))*...
    (1-.5*((X.^2 + Y.^2)/(logsig^2))).*...
    exp((-(X.^2 + Y.^2))/(2*logsig^2));
kern = kern./sum(kern(:));

%% parameters for gaussian filtering:
gausssig = 1;     %in microns
gausssig = gausssig/pix2mic;

%% some pre-processing to smooth out noise
numframe = size(mov,3);
smoothMov = zeros(size(mov));
waitbarHandle = waitbar(0,'filtering movie');
for n = 1:numframe
    temp = mov(:,:,n);
    temp = medfilt2(temp);
    if 0
        %LoG filter; 
        temp = conv2(temp,kern,'same');
        temp(temp<0) = 0;
    else
        temp = imgaussfilt(temp,gausssig);
    end
    smoothMov(:,:,n) = temp;
    waitbar(n/numframe)
end
delete(waitbarHandle);

%% get vessel centers either through automatic detection, a preloaded image, or selection in matlab
if nargin>2 && ~isempty(autoSelectCenters) && autoSelectCenters
    [edgeim, centim] = detectVesselCenters(smoothMov);
    [ep1,ep2,vesselID] = linescansFromCenters([],centim);
    valids = cullEndpoints(ep1,ep2,edgeim);
    ep1 = ep1(valids,:);ep2 = ep2(valids,:);vesselID = vesselID(valids);
    % testing block to make sure the automatic detection is working okay:
    if 1
        figure
        imagesc(smoothMov(:,:,30));
        hold on
        for n = 1:size(ep1,1)
            plot([ep1(n,1),ep2(n,1)],[ep1(n,2),ep2(n,2)],'linewidth',3)
        end
        hold off
    end
elseif nargin > 1 && ~isempty(vesselCenterIm) && ~isstruct(vesselCenterIm)
    % start off with manually selected linescans from FIJI
    [ep1,ep2] = linescansFromCenters(vesselCenterIm);
else
    % alternatively, can select centers directly in Matlab from the image
    if isstruct(vesselCenterIm)
        ep1 = vesselCenterIm.ep1;
        ep2 = vesselCenterIm.ep2;
        vesselID = vesselCenterIm.endpointIDs;
    else
        [ep1,ep2,vesselID] = linescansFromCentersMalabSelect(mov);
    end
end


%% generate kymographs from linescans
waitbarHandle = waitbar(0,'generating kymographs');
kymo = cell(1,size(ep1,1));
for n = 1:length(kymo)
    kymo{n} = linescan2kymo(smoothMov,ep1(n,:),ep2(n,:));
    waitbar(n/length(kymo))
end
delete(waitbarHandle);

%% track kymographs
edgepos = cell(size(kymo));
roughedges = cell(size(kymo));
r = zeros(length(kymo),numframe);
kymfig = figure;
for n = 1:length(kymo)
    try
        [roughedges{n},refinededges] = kymtrack_v3_4(kymo{n},true);     %older iteration used kymtrack_v3_3
    catch
        disp('error on kymograph tracking')
        roughedges{n} = nan(2,size(kymo{n},2));
        refinededges = nan(2,size(kymo{n},2));
    end
    r(n,:) = diff(roughedges{n});
    edgepos{n} = refinededges;
    figure(kymfig)
    imagesc(kymo{n})
    hold on
    plot(edgepos{n}(1,:),'r.');plot(edgepos{n}(2,:),'g.');
    hold off
    pause(.02)
end

temp = mean(r(:,25:45),2);
temp = repmat(temp,1,numframe);
temp = r./temp;
ddd = temp;             %delta d/d sort of thing

%sort through ddd and eliminate bad trajectories
stdwindow = 5;
stdwindowThresh = 2.5;
for n = 1:size(ddd,1)
    cur = ddd(n,:);
    s = movstd(cur,stdwindow);
    s(s>stdwindowThresh) = nan;
    ddd(n,isnan(s)) = nan;
end

results.rawTracking = r;
results.edgepos = edgepos;
results.kymos = kymo;
results.roughedges = roughedges;
results.endpoints = cat(3,ep1,ep2);
results.vesselID = vesselID;